#function for computing the likelihood
function Likelihood(prim::Primitives, est::Estimands, res::Results, estim_data::Array{Float64,2})
    @unpack μ_grid, μs_grid, e_grid, m_grid, p_grid, ℓ_grid, a_grid_1, x_grid_1, ac_grid, f_grid, ℓp_type_grid, h_grid, n_τ, e_s_grid = prim #grids
    @unpack n_μ, n_μs, n_e, n_m, n_p, n_ℓ, n_a_1, n_x_1, n_ac, n_f, n_ℓp, n_h, n_e_s = prim #grid dimensions
    @unpack β, λ_0s_vec, λ_1s_vec, λ_2s_vec, λ_3s_vec, div_chars = prim #discount factor and spouse wage parameters
    @unpack α_1, α_2, α_3, α_4, α_5, α_6, α_7, α_8, τ_s, τ_p, λ_0, λ_1, λ_2, λ_3, γ_0, γ_1, γ_2, γ_3, γ_4, σ_ε, σ_ξ = est
    @unpack γ_5, τ_prob, λ_4, λ_5, σ_θ, λ_6, μ_prob = est
    @unpack emax_1, emax_2, cutoff_1, cutoff_1_fert = res
    
    dist = Normal()
    root_u = sqrt(σ_ε^2 + σ_ξ^2)
    ρ = σ_ε/root_u

    nrow = length(estim_data[:,1]) #store number of rows in data
    likelihood_total = 0.0 #preallocate likelihood value
    lhoods_type = ones(4) #preallocate likelihoods by unobserved het. 2x3 = 4 types total
    #τ_probs = [1-τ_prob, τ_prob]
    #τ_μ_probs = [(1-τ_prob)*0.333, (1-τ_prob)*0.333, (1-τ_prob)*0.333, τ_prob*0.333, τ_prob*0.333, τ_prob*0.333]
    τ_μ_probs = [(1-τ_prob)*(1-μ_prob), (1-τ_prob)*(μ_prob), τ_prob*(1-μ_prob), τ_prob*(μ_prob)]

    #type probabilities for each peson in sample, to be filled in
    nids = length(unique(estim_data[:,1]))
    ntypes = 4
    lhoods = ones(nids, ntypes)
    row_lhood = 0
    id_current = 0

    #begin main loop
    for i = 1:nrow #loop over rows in data
        row = estim_data[i,:] #grab the row
        id = Int64(row[1]) #ID of the person we're on

        #unpack state variables
        i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_f, i_ℓp, i_e_s =
        Int64(row[4]), Int64(row[5]), Int64(row[6]), Int64(row[7]), Int64(row[8]), Int64(row[9]),
        Int64(row[10]), Int64(row[11]), Int64(row[12]), Int64(row[19])

        if i_e_s == 0
            i_e_s = 1
        end

        #choice variables
        h, wage, lprime = Int64(row[13]), row[14], Int64(row[15])

        #unpack some useful stuff
        e, age, ac, m, x, f = e_grid[i_e], i_a + 21, ac_grid[i_ac], m_grid[i_m], x_grid_1[i_x], f_grid[i_f]
        e_s = e_s_grid[i_e_s]
        

        if id!=id_current #onto a new person
            ##########Update likelihood value
            likelihood_total += log(sum(τ_μ_probs.*lhoods_type))
            row_lhood+=1 #move to new row to update
            id_current = id #reset id
            lhoods_type = ones(4) #preallocate likelihoods by unobserved het. 2x3 = 6 types total
        end

        #nab location wage fixed effect
        η_h, δ, coli, η = div_chars[i_ℓp, 2], div_chars[i_ℓp, 5], div_chars[i_ℓp, 4], div_chars[i_ℓp, 9]
        if i_ℓ!=1
            η_h, δ, coli, η = div_chars[i_ℓ-1, 2], div_chars[i_ℓ-1, 5], div_chars[i_ℓ-1, 4], div_chars[i_ℓ-1, 9]
        end
        Δ = γ_0 - γ_1 * e + γ_2 * (age-21) + γ_3 * (ac>-1) + γ_4 * (m>0)  #compute moving cost
        trans_states = [e, m, age, e_s] #things that govern probabilities of marriage, fertility

        i_ac_p = 1 #index of next-period child age state
        if ac>=0 && ac<4 #kid aged  0 to 3
            i_ac_p = i_ac+1 #next period: kid aged one year older
        end
        if f == 1 #check fertility status
            i_ac_p = 2 #kdi aged 0 next period if woman is pregnant
        end
        lhood_column = 0 #column in likelihood to update: 1-1, 1-2, 2-1, 2-2

        for i_τ = 1:2, i_μ = 1:n_μ #loop over unobserved grandparent type type
            lhood_column+=1

            #unobserved wage type
            μ = μ_grid[i_μ]
            wage_predict = λ_0 + λ_1*e + λ_2*x + λ_3*x^2 + μ + η * (1 + λ_6*e) + λ_4 * (ac==0 || ac==1) + λ_5 * (ac>1)  #agent wages

            #cutoffs
            cutoff = cutoff_1[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_f, i_ℓp, i_τ, i_e_s]
            cutoff_fert = cutoff_1_fert[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s]

            #pritnln(cutoff, " || ", cutoff_fert)

            prob_f = cdf(dist, cutoff_fert/σ_θ)
            if f == 1
                prob_f = 1 - cdf(dist, cutoff_fert/σ_θ)
            end

            #probability of working
            prob_h = cdf(dist, cutoff/σ_ε) #assuming not wokring
            if h == 1 #update lfp; working. Add in measurement error stuff

                #step 1: compute total error
                u = log(wage) - wage_predict

                #assemlbe probability
                frac = (cutoff-(ρ^2)*u)/(σ_ε*sqrt(1-ρ^2)) #assembly of probability
                prob_1 = 1-cdf(dist, frac) #assembly of probability
                prob_2 = (1/root_u)*pdf(dist, (u/root_u)) #assembly of probability
                prob_h = prob_1*prob_2 #formation of probability
            end

            #probabilty of location decision
            i_p_p = h + 1 #next-period past-participation state
            i_x_p = min(i_x + h, n_x_1)
            i_x_p_2 = i_x + h

            prob_ℓ = 0
            if lprime!=99 #ignore if missing
                if i_a == n_a_1 #in terminal age of phase 1
                    emax_grid = zeros(n_ℓ+1)

                    #parent location option
                    emax_grid[1] = β * emax_2[i_μ, i_e, i_m, i_p_p, 1, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)

                    #moving options
                    for i_ℓc = 2:(n_ℓ)
                        emax_grid[i_ℓc] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓc, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                    end

                    #stay8ing option
                    emax_grid[n_ℓ+1] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓ, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s]

                    prob_ℓ = exp(emax_grid[lprime]) / sum(exp.(emax_grid))
                    if lprime == 11 && i_ℓ == 1  #double likelihood if chose to stay in parent location
                        prob_ℓ += exp(emax_grid[1]) / sum(exp.(emax_grid)) #add on probability of parent option, since we can't distinguish the two
                    end
                elseif i_a!=n_a_1 #not in terminal age
                    emax_grid = zeros(n_ℓ+1)

                    #loop over potential marriage/fertility realizations
                    for i_mp = 1:n_m, i_e_s_p = 1:n_e_s
                        mp = m_grid[i_mp] #standardize
                        e_s_p = e_s_grid[i_e_s_p]
                        prob = Compute_Prob_Marr(mp, e_s_p, trans_states, prim) #compute probability of realization. Slight bottleneck here.
                        emax_grid[1] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, 1, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)) * prob
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[i_ℓc] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, i_ℓc, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])) * prob
                        end
                        emax_grid[n_ℓ+1] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, i_ℓ, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p]) * prob
                    end

                    prob_ℓ = exp(emax_grid[lprime]) / sum(exp.(emax_grid))
                    if lprime == 11 && i_ℓ == 1  #double likelihood if chose to stay in parent location
                        prob_ℓ += exp(emax_grid[1]) / sum(exp.(emax_grid)) #add on probability of parent option, since we can't distinguish the two
                    end
                end
            end

            if lprime == 99
                prob_ℓ = 1 #don't use moving likelihood if unobserved
            end

            if age>=25 || e == 0 #check for age
                #println(prob_ℓ, " || ", prob_h, " || ", prob_f)
                lhoods_type[lhood_column] *= (prob_ℓ * prob_h * prob_f) #update the type-spedcific likelihood for the person we're looking at
                lhoods[row_lhood, lhood_column] *= (prob_ℓ * prob_h * prob_f) #update individual likelihood
            end
        end

    end
    likelihood_total += log(sum(τ_μ_probs.*lhoods_type)) #add on the likelihood of the very last person
    likelihood_total, lhoods #return likelihood
end

#objective function to minimize
function Objective_function(g::Array{Float64,1}, dir::String; race::Int64 = 1)
    obj = 0.0 #preallocation of objective function value

    println("Current Guess")
    println(round.(g, digits = 4))

    #load appropriate utilities
    package = Readin(dir)    
    if race == 2
        package = Readin_whites(dir)
    elseif race == 3
        package = Readin_blacks(dir)
    end
    estim_data, trans_probs, div_chars = package[1], package[3], package[4]

    if g[1]<0 || g[2]<0 || g[3]>0 || g[5]<0 || g[6]<0 || g[10]<0 || g[11]<0 || g[11]>1 || g[17]<0 || g[18]<0 || g[12]<0 || g[12]>1 || g[25]>1 || g[25]<0 || g[35]<0 || g[47]<0 || g[49]>1 || g[49]<0
        obj = Inf
    else
        prim, est, res = Initialize(g, trans_probs, div_chars) #initialize primitives, estimands, and results vectors
        Backward_Induct(prim, est, res; race = race) #compute all value functions
        likelihood, lhoods = Likelihood(prim, est, res, estim_data)
        obj = -1 * likelihood

        #garbage collection        
        finalize(res.emax_1)
        finalize(res.cutoff_1)
        finalize(res.emax_2)
        finalize(res.cutoff_2)
        finalize(res.emax_3)
        finalize(res.cutoff_3)
        @everywhere GC.gc()
    end

    #report how we're doing
    println("")
    println("Current error: ", round(obj, digits = 4))
    println("")
    obj #return deliverable
end
###############
